

		clear
		use Data_20170719.dta
		set seed 5500

			/* 
			Note that the marginal effect is set to increase from 500 to 1500 troops.
			This represents a change from 500 below the median non-zero value of deployments (~1000 troops) 
			to 500 above the median. 500 is also about the size of one US battalion. 
			*/
			
		* Marginal Effect of an increase in peacekeepers across duration of US presence.
			levelsof duration_us, local(levels)	
			xtset ccode year
			probit pdest c.duration_us##c.ln_troops intensity_high w_intensity_high w_psource w_pdest polity2 ln_pop ln_gdpwb  border, robust 

			sum duration_us if e(sample)
			local min = r(min)
			local max = r(max)
			local interval = ((`max'-`min')/50)
			margins, at(ln_troops=6.22 ln_troops=7.32 duration_us=(`min'(1)`max') (mean) w_psource w_pdest polity2 ln_pop ln_gdpwb border w_intensity_high intensity_high=0) mcompare(bonferroni) vsquish level(90) post  		

			local a = 4
			local b = 3
			#delimit ;
			matrix list e(b);
			mlincom 2-1, clear stats(all)	;
			while `a'<=60 { ;
			mlincom `a'-`b', stats(all) add ;
			local a=`a'+2 ;
			local b = `b'+2 ;
			} ;
			matrix list _mlincom			;
			svmat double _mlincom, names(col);
			range xvar `min' `max' 15 ;
			gen yline = 0 ;
			
			#delimit ;
			twoway rarea ll ul xvar,  color(gs12%50) lcolor(gs12%0) fi(100)
			|| line lincom xvar, lpattern(solid) msym(O) lwidth(.3) 
			yline(0)
			xlabel(0(5)65, format(%9.0fc)) 
			ylabel(, format(%9.2fc))
			xtitle("Duration of US Troop Presence")
			ytitle("Change in Predicted Probability")
			legend(off)
			title("")
			aspect(1)
			saving(me_us_duration, replace)
			;
			#delimit cr
			
			graph export "me_pdest_us_duration.pdf", replace
		
			
